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We present an alternative to the conventional matrix product state representation, which allows 
us to avoid the local Hilbert space truncation many numerical methods employ, as well as drastically 
reduce the number of matrices required to describe the state. Utilising chain mappings (linear as 
well as logarithmic) of the spin-boson model (SBM) Hamiltonian onto a semi-infinite chain, we apply 
the new method to the sub-ohmic SBM, where we can reproduce many well established features of 
the quantum phase transition, such as the critical exponent | predicted by mean-field theory. Via 
extrapolation of finite-chain results we are able to determine the infinite-chain critical couplings a c 
at which the transition occurs. 



INTRODUCTION 

The spin-boson model (SBM) describes a single two- 
level system (TLS), a spin, coupled to environmental 
degrees of freedom represented by a continuos bath of 
bosonic field modes. It is one of the most important 
models for studying the general effects arising when a 
quantum system is coupled to an environment [I] . In the 
sub-ohmic version, it possesses a mean-field like quantum 
phase transition between a localised and a delocalised 
phase at zero temperature. This phase transition has 
been the subject of extensive numerical and analytical in- 
vestigations [2HE]. Yet, many numerical approaches face 
challenges near and above the transition to the localised 
phase. This is due to the rapidly rising number of field 
excitations in the localised phase, which imply that the 
quantum states of the field modes span an increasingly 
large subspace of their full Hilbert space. Most numer- 
ical methods are however based on local Hilbert space 
truncation and are hence discarding more and more vital 
information. 

In this paper, we describe a variation on the matrix 
product state (MPS) representation that avoids the local 
Hilbert space truncation. We demonstrate its viability 
by applying it to the sub-ohmic SBM to study the prop- 
erties of the quantum phase transition and compare them 
to an analytical approach based on a variational ansatz 
0- 



MPS METHOD 

Generally, the compound quantum state of an M-level 
system (labelled by k = 1, ...,M) and an environment 
of N bosonic modes (labelled by i\, ...,ijy G N°) can be 
described by 

\ip) = ^2 c k ,i u ...,i N \k,ii,...,i N ) . (1) 

k,i± : ...,iN 



Alternatively, it is possible to express the coefficients 
c~k,i u ...,i N as a product of \ x \ matrices, where the sys- 
tem as well as each mode have a unique set of matrices 
associated with them 

TV N 

k,ii,...,iN m— 1 rn — 1 

N 

= E tT l Sik) II B<My/iJ]\k,iu...,iN), (2) 
fe,il,...,ijv m— 1 

where |0) represents the vacuum state of all modes and 
&Tj is the creation operator of the m-th mode. 
This representation of the state is known as a matrix 
product state (MPS) [HHH3- In this form the state of 
the M-level system is represented by the M matrices 
S^, S( M "> (e.g. M=2 for a qubit, like in the SBM) 
and the m-th field mode is represented by a semi-infinite 
set of matrices {Bm" 1 '}. In these infinities lies the prob- 
lem of the traditional MPS. Numerical calculations are 
limited to a finite set of matrices, hence the local Hilbert 
space associated with each mode has to be truncated to 
a finite size by limiting the local dimension. Under cer- 
tain conditions such as high mean excitation numbers, 
this truncation can lead to substantial errors in numeri- 
cal calculations. 

Here we propose an alternative MPS that avoids trunca- 
tion. This is achieved by reducing the number of matrices 
per mode to a single matrix X m , which is defined such 
that 

B (i m ) = (3) 

where < i m < oo, i.e. the infinite set {Bm m ^} is now 
formed from powers of a single matrix X m , reducing the 
total number of matrices required to fully describe the 
state to M + TV [13]. The additional factor of (i m !) _1 is 
chosen to simplify later calculations. 
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Thus, substituting into ([2]), the MPS can be written as 

X 



N 



e e M^n^i^^- 



an) 



(4) 



k— 1 i\ 



...i N =0 



This form not only enables us to avoid the truncation, 
but also makes it quite straightforward to determine the 
normalization as well as expectation values. Introducing 
the new notation 



-.6 _ 



n 



,X m (g>X„ 



and 



M 



k=l 



we find for the norm 



(0#) = tr[ET* 



Ni 



(5) 



(6) 



(7) 



The exponentiation of the X m matrices in combination 
with deliberate choice of the (im!)^ 1 factor results in the 
exponential functions, which are straightforward to eval- 
uate numerically. As an example of an expectation value, 
we find for the population of the fc-th mode 



(VlWV>) = tr[ET k 1 (X k ® X k )T 



fc+u 



(8) 



In general we find that all quantities of interest are simple 
traces over products of x 2 x x 2 matrices. 
The MPS allows us to determine the ground state of a 
system described by an arbitrary Hamiltonian H as the 
state \tp) for which the S and X minimize the energy 



E = 



(9) 



We now apply this method to analyze the ground state 
properties of the sub-ohmic spin-boson model. 



SPIN BOSON MODEL 

The Hamiltonian of the (unbiased) spin-boson model 
(SBM) is given by (h = 1) 

H S b = -\&°x + ^<?z y^9i(ai + a\) + E w l a \ a l> ( 10 ) 



where Oi are the usual Pauli matrices describing a two- 
level system (TLS) with tunnelling amplitude A. ai and 
at are the bosonic annihilation and creation operators 
of the environment, which consists of bath modes with 
frequency u>i. The key quantity in the description of the 
system-environment interaction is the spectral function 



J(ui) — T^^2ig 2 5(uj — u>i). Here we consider a spectral 



function of the form 



(11) 



as given in |14) . where uj c is the maximum cut-off fre- 
quency of the spectrum. In the following we focus ex- 
clusively on the sub-ohmic case for which < s < 1, in 
particular on the case s < 0.5. 

In [21 [TS] it was shown that a Hamiltonian of the form 
eq. ( 10 ) can be mapped onto a semi-infinite chain of 



bosonic modes that experience nearest neighbour inter- 
action only, with the system only coupling to the first 
chain site. The transformed Hamiltonian can be written 



as 



H = - -Acr^ 
2 



c cr z (6o + &o) 



+ V WrJ&bn + t n (bl +1 b n + bib n +i), (12) 



n=l 



where the coupling strength between the TLS and the 
first site is given by 



Co 



2(s + 1) 



(13) 



and the local energies and tunnelling amplitudes of the 
sites are 



1 



(s + 2n)(2 + s + 2n) 



and 



(jj c (l + n)(l + s + n) 3 + s + 2n 



s + 2 + 2n) (3 + s + 2n) V 1 + s + 2n 



(14) 



(15) 



respectively. This mapping brings several advantages 
which include the analytical forms for the parameters of 
the resulting chain model given above, an intuitive pic- 
ture how irreversibility emerges |16j and the ready appli- 
cability of the MPS method 

Utilizing the MPS eq. Q, we can now find the ground 
state of the spin-boson system by minimizing eq. d9]). 
Specifically, we have to find the S and X which minimize 
(assuming ('01'!/') = 1) 



E = (ip\H\<iP) = E ioc + E c + E, 



chain j 



(16) 



with 

Eloc 



E c = c (a z (b + bl)), 

OO 

Echain = W n (b^bn) 



(17) 
(18) 

t n ((bi +1 b n ) + (btb n+1 )). (19) 



In terms of the MPS description we get after truncating 
the chain length to N sites 



3 



E 



2 

JV-l 



5^ {ir [sTJ (w„(X„ ® X„) + i„(A„ +1 <g> X n + X n ® X n+1 ))r^ +1 ] } + Wjv tr [Hrf (Xjv ® a w ) 



(20) 



This minimization can be done numerically to yield the 
full ground state MPS of the TLS and the (truncated) 
chain. In the following we will present some results for 
the ground state properties of the sub-ohmic SBM. The 
minimizations were carried out using 2x2 matrices and 
MATLAB's fminunc function. 



RESULTS 

The sub-ohmic SBM possesses a mean-field like 
continuous phase transition in the magnetisation for 
< s < 0.5 at a critical coupling strength a c between 
system and environment. For small coupling strengths 
a < a c the TLS is in a delocalised phase, having no net 
magnetisation. Above the critical point a > a c the en- 
vironment induces a spontaneous magnetisation on the 
TLS, which then exhibits a doubly degenerate localised 
phase. 

In the delocalised phase the mean site population is ex- 
pected to diverge along the chain. This has so far made it 
difficult to accurately study the system close to and above 
the phase transition with numerical methods that rely on 
Hilbert space truncation, such as numerical renormalisa- 



0.20 




FIG. 1. Critical coupling a c as a function of the inverse chain 
length 1/N for A = 1 and ui c = 10. The fitted functions 
are of the form a c = ae b ^ N , where a is the extrapolated lim- 
iting value for a c as N — > oo. The dotted line for s — 0.5 
corresponds to logarithmic discretisation (A = 1.5) and the 
short dashed lines are the critical values a c predicted by the 
polaron ansatz in [2J. 



tion group (NRG) [3] and density matrix renormalisation 
group (DMRG) methods [T7J [T5] . We propose that this 
truncation and many of the associated problems can be 
avoided using the MPS of form eq. Q and the minimiza- 
tion eq. poL 



However, one still has to perform a different kind of 
truncation to make the numerical simulation feasible, 
namely truncation of the chain length N which effec- 
tively amounts to neglecting low frequency components 
of the environmental bath. 

Fig. [T] shows the critical coupling strengths a c plotted 
against the inverse chain length 1/N for several values 
of s, which we determined using the new MPS method 
described above (with 2x2 matrices). A good fit to the 
data was found for an ansatz of the form a c — ae b ^ N , 
which we then fitted for each s and were thus able to ex- 
tract an extrapolated value for a c in the limit N — > oo. 
Table [T] shows the results of the extrapolation. In [5] a 
variational ansatz was used to predict some of the proper- 
ties of the sub-ohmic SBM ground state. For the critical 
coupling strength they predict a value 



sin(7rs)e s / 2 / A ^ 1 ~ 
27r(l — s) \uj c 



(21) 



These predicted values are indicated in Fig. [T]by dashed 
lines and are also listed in Table |T] along with the frac- 
tional deviation between our extrapolation results and 
the predicted values. We find that the predictions agree 
well with our result for large s but show significant de- 
viations at smaller s. This suggests that the mean field 
type ansatz for the environment holds well for larger s 
but fails for decreasing s, possibly as a result of the in- 
creasing correlations in the environment. 
Being able to find an MPS representation for the ground 



TABLE I. Critical coupling strengths a c for various values of 
s, together with the respective critical values a c derived from 
the polaron ansatz in [2] and the fractional deviation. 



s 


Ctc 




(a c — a c )/a c 


0.1 


0.0132 


0.0065 


1.014 


0.2 


0.0216 


0.0168 


0.286 


0.3 


0.0353 


0.0316 


0.117 


0.4 


0.0552 


0.0519 


0.064 


0.5 


0.0828 


0.0784 


0.056 
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state, we were also able to analyse the general proper- 
ties of the state in both the delocalised and the localised 
phase. Fig. [2] shows the magnetisation M — \ (a z ) | of 
the TLS for some representative values of s. This and 
the following results were obtained with a chain length 
N = 50. The figure clearly shows the two phases, sep- 
arated by a second-order transition. In the delocalised 
phase a < a c the order parameter M is zero. Above 
the critical coupling strength a c the TLS obtains a finite 
magnetisation with a tendency to full localisation M = 1 
as a grows large. This localised phase is two-fold degen- 
erate with M — ± (a z ) both being solutions. 
Mean-field theory predicts a second-order magnetic tran- 
sition at a c with 



Moc \a c -a\ 1/2 . 



(22) 



Our simulations do indeed reproduce the correct critical 
mean-field exponent 1/2 well, as can be seen in Fig. [3j 
where we have plotted the magnetisation for a > a c on 
a log- log-plot. The mean-field result is indicated in the 
figure by the solid straight line. 

In |] a variational ansatz was used to predict the 
amount of entanglement in the TLS, defined as the von 
Neumann entropy of its reduced density matrix. Our 
results are presented in Fig. [4] Despite the slight de- 
viations in the values for a c , which most likely arise 
as a combination of the finite chain length we consid- 
ered as well as the inherent differences between the two 
approaches (c.f. Table ttj), our results are in excellent 
qualitative agreement with their predictions. In the de- 
localised phase entanglement increases. At the critical 
coupling a c the entanglement exhibits a cusp and then 
decays rather rapidly in the localised phase due to the 
system evolving into a product state. 
In addition to the entanglement of the TLS, we also 
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FIG. 2. Magnetisation M as a function of a for N = 50, 
A = 1 and ui c — 10. The first order discontinuity in the order 
parameter M at a c marks the phase transition. The magenta 
points in this and the following figures correspond to s = 0.5 
with logarithmic chain discretisation. Data points are joined 
for better visibilty. Whereas in [2] similar plots could only 
be found through an analytical ansatz, the new method now 
allows us to obtain numerical results. 
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FIG. 3. Log-Log plot of Magnetisation M as a function of 
(a — a c )/a c for a > a c , N = 50, A = 1 and oj c = 10. The 
solid straight line represents the expected mean field exponent 
1/2 just above the critical coupling a c . 



looked at the entanglement of the individual sites in 
the chain. As a representative example, the results for 
s = 0.3 are shown in Fig. [5] for the first ten chain 
sites. We find that only the first few sites carry sig- 
nificant amounts of entanglement and that the general 
behaviour with changing a closely resembles the entan- 
glement properties of the TLS in Fig. |4j A substantial 
spread of entanglement along the chain is only observable 
very close to the phase transition. 

Another observable of interest is the coherence (<r x ) 
which is shown in Fig. [6] as a function of a. The re- 
sults are again in excellent qualitative agreement with 
the results from the variational approach [2]. The co- 
herence is continuously decreasing, with a faster decay 
above the transition, at which we observe a cusp. 
As mentioned above, the main reason why other numeri- 
cal approaches failed to return accurate results near and 
above the critical coupling is that, whereas in the delo- 
calised phase the mean occupation {b' n b n ) of chain site 
n rapidly decreases along the chain, it considerably rises 
in the localised phase. In fact, 2j predicts that (b^bn) 
diverges along the chain for a > a c . In Fig. [7] we plot 
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FIG. 4. Entanglement (von Neumann entropy) between the 
TLS and the environment for N = 50, A = 1 and oj c = 10. 
The maxima coincide with the phase transition at a c . 
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FIG. 5. Entanglement E of the individual chain sites as a 
function of a for the first 10 sites for s = 0.3. The entangle- 
ment of site n is here defined as the von Neumann entropy 
of the reduced density matrix of site n. Only the first few 
sites show significant amounts of entanglement. This shows 
that the analytical ansatz in [2] does have a sound basis in 
this regime but fails to be accurate near the critical point and 
perhaps for other choices of s when entanglement is larger. 



(b\b n ) as a function of a and n. Below the transition 
we find that indeed the average population of the sites 
decreases with n. At the transition we observe a sudden 
increase in (b^bn) and the maximum begins to shift away 
from the first site, further along the chain. This rapid rise 
in the occupancy shows why methods such as DMRG, 
which rely on Hilbert space truncation, are challenged 
in this regime, since the information about the system 
is spread over more and more basis states, only a finite 
number of which these methods retain. 
An alternative method to the linear chain mapping we 
have thus far considered, is provided by logarithmic dis- 
cretisation of the spectrum [31 [U [TH [TO] , which is exten- 
sively used in NRG. It does not linearly subdivide the 
bath spectral function, but instead splits it in intervals 
[A~(™ +1 ), A - "], where A > 1 is the discretisation param- 
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FIG. 6. Expectation value (cr x ) as a function of a for N — 50, 
A = 1 and uj c — 10. A cusp at a c marks the phase transition. 



FIG. 7. Expectation value (6„6n) as a function of a for each 
chain site (labelled by n) for s — 0.3. Above the phase transi- 
tion the mean population of the sites quickly grows, rendering 
numerical methods such as DMRG which are based on Hilbert 
space truncation inaccurate in this regime. 



eter and n <G N°. This new Hamiltonian has again the 
same form eq. (12), but with different site frequencies 



and transition amplitudes 



C n ) 



and 



_ Nn ± l A 



(23) 



(24) 



respectively, where Q s , A n , C n and N n are given in [3j. 
To be applicable for numerical methods it is again nec- 
essary to truncate the resulting chain Hamiltonian at a 
finite number of sites n = N. 

A challenge for many numerical methods with the loga- 
rithmic discretisation is the fact that the mean occupa- 
tion of the chain sites is on average considerably larger 
than on the linearly discretised chain. This quickly leads 
to the breakdown of these methods. However, with our 
new method avoiding the truncation, we were again able 
to establish numerical results for the sub-ohmic SBM 
ground state. The dotted line in Fig. [I] shows the same 
extrapolation for s = 0.5 as done before, now using the 
logarithmically discretised Hamiltonian with discretisa- 
tion parameter A = 1.5. Despite the missing data for 
large N, we see that the critical coupling converges to 
a similar value as found before, with a fractional devia- 
tion (a c - a [ c LD) )/a[ LD) = 0.030. We also find that with 
the same chain length N, the logarithmically discretised 
Hamiltonian results in a value for a c that is closer to the 
limiting value for N — > oo found via extrapolation than in 
the linearly discretised case. However, for larger values 
of N, the simulation takes considerably longer to con- 
verge near the phase transition, hence we did not acquire 
reliable data points for N > 25. The dotted magenta 
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lines for s = 0.5 in Figures [2j [3j [4] and [6] also show re- 
sults using the logarithmic discretisation with N = 50, 
using exactly the same computational time as the results 
for linear discretisation. The results are almost identical 
to those obtained via linear discretisation, except just 
above the transition at around 0.09 < a < 0.11 where 
the convergence issues come into play. 
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CONCLUSION 



By modifying the traditional matrix product state rep- 
resentation, we were able to avoid the local Hilbert space 
truncation that leads to the failure of many numerical 
methods in a regime of high field mode excitation. Us- 
ing this modified state representation combined with a 
method of energy minimisation, we were able to give a 
detailed study of the ground state properties of the sub- 
ohmic SBM. Our findings are in good agreement with 
previous numerical and analytical results, but extend 
these to new regimes of the spin-boson model, partic- 
ularly the region close to and above the phase transi- 
tion, which has so far not been accessible to numerical 
methods, since methods such as DMRG fail to produce 
reliable results due to the rapid increase of the Hilbert 
space dimension of the environmental modes. In addi- 
tion, our method allowed us to give an analysis of the 
chain properties such as mode excitation and entangle- 
ment for specific sites along the chain. It also has the 
advantage of being comparatively easy to implement nu- 
merically. Hence the method provides a promising new 
tool to investigate the localised phase of the SBM near 
and above the transition and to test the current analyti- 
cal results such as the mean field type approaches in this 
regime. 
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